Análisis de Series de Tiempo de Homicidios en colombia

Chapter 1 Descripción

La dinámica y evolución de los homicidios en Colombia constituye un fenómeno de gran relevancia que merece ser minuciosamente analizado y comprendido. Pues, el entendimiento de sus tendencias, patrones y factores resulta esencial para abordar la complejidad de este delito en el país. En este sentido, este análisis podría arrojar luz a aspectos cruciales como la variación estacional de los homicidios, los potenciales impactos de factores socioeconómicos y demográficos y/o la influencia de las políticas de seguridad implementadas en el territorio nacional.

La base de datos proporcionada tiene un espacio temporal que data desde el 2010 hasta enero de 2024.

Fuente: Dirección de Investigación Criminal e Interpol (DIJIN) - Policía Nacional de Colombia.

GOV.CO_DatosAbiertos

Análisis exploratorio

Dimensiones

Se realiza la exploracion de las dimensiones de la base de datos evidenciando que esta cuenta con:

Filas: 75152 Columnas: 9

Tipo de variables

Se debe corrigir el tipo de las columnas FECHA HECHO y CANTIDAD, dado a que éstas son de tipo Date y Número, luego de aplicar los cambios se observa:

CANTIDAD: La variable es tipo numeric
FECHA HECHO: La variable es tipo Date

Identificación de registros vacíos

Se realiza la verificacion de valores nulos o vacios obteniendo los siguientes resultados:

##         DEPARTAMENTO            MUNICIPIO          CODIGO DANE 
##                    0                    0                    0 
##         ARMAS MEDIOS          FECHA HECHO               GENERO 
##                    0                    0                    0 
##         GRUPO ETARÍO DESCRIPCIÓN CONDUCTA             CANTIDAD 
##                    0                    0                    0

Clara mente se observa que no existen valos nulos o vacios en ninguna de las variables.

Identificación de registros vacíos

Se realiza la verificacion de valores nulos o vacios obteniendo los siguientes resultados:

##  /\     /\
## {  `---'  }
## {  O   O  }
## ==>  V <==  No need for mice. This data set is completely observed.
##  \  \|/  /
##   `-----'

##       DEPARTAMENTO MUNICIPIO CODIGO DANE ARMAS MEDIOS FECHA HECHO GENERO
## 75152            1         1           1            1           1      1
##                  0         0           0            0           0      0
##       GRUPO ETARÍO DESCRIPCIÓN CONDUCTA CANTIDAD  
## 75152            1                    1        1 0
##                  0                    0        0 0

Clara mente se observa que no existen valores nulos o vacios en ninguna de las variables.

- Resumen de Estadisticos

A continuacion se muestra un resumen de los estadisticos basicos de nuestra variable CANTIDAD:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   1.000   1.000   1.105   1.000  16.000

- Transformación del dataset

A continuacion presentamos un resumen de los estadisticos de los homicidios totales de la serie de tiempo luego de realizar la transformacion de los datos:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   174.0   397.0   449.0   491.5   526.0  1264.0

Dado a que es una serie de tiempo solo tendremos en cuenta la construcción de una base que contenga las variables FECHA HECHO y CANTIDAD para el analisis a realizar, de igualmanera se determinan los siguientes parametros:

Frecuencia de la serie: Anual
Inicio de la Serie: 2010, 1
Fin de la Serie: 2024, 1

adicionalmente se resaliza la verificacion de la clase de la serie de tiempo:

Calse de la serie: ts

- Garficos

Luego de analizar los resultados se evidencia que en el mes de abril hay mayor numero de homicidios, adicionalmente se evidencia que la media de los meses se encuentra entre los 400 y 500 homicidios.

En 7 meses se observan unos valores atipicos que superan los 1000 homicidios y uno en el que se presentaron menos de 200, seria de gran valor hacer un analisis detallado de estos datos con el objetivo de entender mejor la naturaleza de estos resultados.

Para el caso de la grafica de rezagos se puede afirmar que no existe aleatoriedad, debido a que no se reflejan patrones identificables en los datos.

- Media Movil

A continuacion se realiza el calculo de las medias moviles (SMA y EMA) de la serie de datos con el objetivo de obtener de forma mas clara el comportamiento de nuestra serie.

- Gráfica de medias moviles exponencial vs simple

Durante los últimos 13 años, los homicidios en Colombia han experimentado un aumento gradual. Las medias móviles de 12 meses muestran que en 2010 había entre 230 y 240 asesinatos, comparados con 1000 a 850 asesinatos en los últimos meses de 2023 y enero de 2024, quintuplicando así las cifras de este fenómeno en el país. Se observa una tendencia a la baja al finalizar el primer semestre de cada año, seguido por un aumento durante los últimos meses, adicionalmente, se identifican dos períodos de fluctuaciones significativas:

Una baja notable al comienzo de la pandemia en 2020, dada la crisis sanitaria provocada por el COVID y la política de aislamiento social
Un aumento sostenido en casi todo 2023, este comportamiento podria estar asociados a aumentos de bandas criminales y grupos armados como efecto de los cambios politicos que se generaron con el actual gobierno en materia de seguridad.

En cuanto a las líneas móviles exponenciales versus las simples, aunque no coinciden exactamente en su posición, sí lo hacen en cuanto a su tendencia, siendo la línea simple más suavizada que la exponencial.

En congruencia con las medias móviles, se observa que la cantidad de homicidios no supera los 375 casos mensuales mensuales antes del 2020, sin embargo en el los periodos posteriores como en el 2023 se observa que se alcanzan valores tope hasta de 1000 muertes mesuales en el país a causa de los homicidios.

- Transformaciones básicas Series de Tiempo

A continuacion realizaremos algunas transformaciones que nos permitiran detallar mucho mejor el analisis de nuestra serie de tiempo.

Acorde a la tendencia el comportamiento de los homicidios es lineal durante los años 2010 y 2020. Posteriormente, se vuelve creciente hasta finales del 2023.

- Identificación de Estacionaridad Dicker - Fuller

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Data_ts$homicidios
## Dickey-Fuller = -4.1312, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Dado que el p-value es menor al nivel de significancia de 0.05 se acepta la hipotesis alternativa de que la serie sí es estacionaria

- Estacionalidad por año

Se evidencian picos en la mayoría de los años principalmente en los meses de febrero, abril, junio y octubre, a excepción de 2023; a diferencia de los meses mayo, septiembre y noviembre, don de la cantidad de hpmicidios disminuye.

- Diferenciación

Dado a que en los modelos de series de tiempo se requiere tener en cuenta la estacionariedad, para una mejor modelización y capacidad predictiva se procede a obtener las diferencias para hallarla. En otras palabras, se realiza para la estructuración del modelo a realizar.

¿Cuántas diferencias se necesitan para hallar estacionariedad?

Despues de realizar el procedimiento de diferenciacion se llega a la conclusion que se solo se requieren 1 diferencias para identificar la estacionariedad.

- Transformación para la variabilidad

A continuacion, se aplica una transformacion logaritmica a la serie de tiempo, esto se realiza para cumplir con el supuesto de que la serie tiene variabilidad constante, para una mayor estabilidad e interpretación de datos.

- Aplicación de diferenciación

- Aplicación de función ACF

El realizar la autocorrelacion nos permite identificar un comportamiento estacionario con respecto al tiempo en la serie de tiempo.

- Aplicación de función PACF

Al identificar la estructura autorregresiva en la serie, se tiene que, se necesitan 3 rezagos para predecir el valor actual de la serie.

Aplicación de Holt-Winters Model.

Dada la evidencia de existencia de un valor medio, tendencia y estacionalidad en los datos; se permite la aplicación del modelo Holt-Winters; como modelo predictorio de largo y mediano plazo por medio de un triple suavizado exponencial al tener en cuenta los aspectos mencionados con anterioridad.

La primera decisión radica en elegir el tipo de patrón de estacionalidad, es decir, si este modelo debería ser representado por una estacionalidad aditiva o multiplicativa. Basándonos en la evidencia previamente encontrada, se observa una tendencia que aumenta o disminuye proporcionalmente con el promedio móvil de los datos, lo que sugiere que no permanece constante. Además, la amplitud de esta estacionalidad varía con el nivel encontrado. Por lo tanto, se decidió que el modelo a escoger será multiplicativo.

Multiplicative

La gráfica nos brinda una visualización del comportamiento de las predicciones generadas por el modelo multiplicativo Holt-Winters. En un primer vistazo, podemos observar un ajuste cercano entre las predicciones y los datos base. Además, se aprecia una tendencia y estacionalidad similar entre ambas.

Los valores de los parámetro obtenidos son: Alpha: 0.694986 Beta: 0 Gamma: 0.0054929

Al revisar los parámetros usados se puede apreciar que el modelo está dando más peso al parámetro Alpha que corresponde a la tendencia.

Para un mejor ajuste, se probará modificar los parametros Beta, Gamma y Alpha del modelo holt-winters, dado a que éstas se asocian respectivamente con tendencia, estacionalidad y nivel (promedio móvil).

Alpha es un número entre 0 y 1 que determina cuanto peso se le da a las observaciones más recientes al calcular la tendencia.

Beta es un número entre 0 y 1 que determina cuanto peso se le da a las observaciones más recientes en cuanto su estacionalidad.

Gamma es un número entre 0 y 1 que determina cuanto peso se le da a las observaciones más recientes en cuanto su nivel (promedio móvil).

Al modificar dando más peso al parámetro Gamma correspondiente al Nivel, se puede apreciar como las estimaciones cambian en magnitud y los picos y valles no coinciden en igual proporción como en el modelo original.

Al modificar dando más peso al parámetro Beta correspondiente a la Estacionalidad, se puede apreciar como las estimaciones del 2024 es muy diferente a los modelos anteriores, lo cual se explica al considerar lo ocurrido en el año 2023 como consecuencia de factores externos que sólo afectan en ese periodo de tiempo.

El ajuste del modelo muestra una tendencia constante, un nivel que se aproxima a los datos observados y una estacionalidad periódica entre los años.

En los valores predichos bajo el modelo Holt-winters, se implementó la predicción de 11 meses posteriores a enero 2024 con intervalos de confianza del 80 y 95%, observando unas barreras correspondientes a 1200 maximo y menos de 200 homicidios sin 0, por observación, el número estaría entre 180 a 150 homicidios.

- Forcast evaluation

Realizando una autocorrelación a través de la función ‘acf’, se establecen unos rezagos de 30, haciendo alusión a la misma obsevación en el mes anterior, y dejando denotado no tener en cuenta los NA’s en la operación.

Dado lo anterior, se observa que solo un error sobresale de la franja inferior de confianza, denotando la existencia de factores externos que tienen un grado elevado de significancia en la variabilidad de los datos pudiendo atribuirse a los cambios políticos del nuevo gobierno central o la crisis sanitaria.

- Testeo de errores

La prueba Ljung-box evalua la hipotesis nula de que no hay autocorrelación en los datos hasta el rezago especificado, en este caso 30.

## 
##  Box-Ljung test
## 
## data:  Modelo_evaluation$residuals
## X-squared = 37.565, df = 30, p-value = 0.1613

Dado que p-value es mayor a 0.05 se acepta la hipotesis nula advertida en el parrafo anterior, lo que significa que en primera instacia hay una buena especificación del modelo, las predicciones son confiables y los residuos son aleatorios, lo que permite al modelo capturar adecuadamente la estructura de dependencia temporal de los datos.

Así las cosas, graficando los residuos del modelo, se tiene que gran parte de las observaciones se encuentran alrededor de 0, lo que sugiere que el modelo no tiene un sesgo sistemático en sus predicciones; captura bien la estacionalidad y la tendencia de los datos; y explica en mayor parte la variabilidad de los datos.

- Modelaje Box - Jenkins

La prueba dicker-fuller realizada en líneas de código anteriores, confirma estacionariedad en los datos con un p-value de 0.01; se ajustó la variabilidad y se halló los rezagos correspondientes. De esta manera, con todas estas observaciones realizadas, este subcapitulo de modelaje se trabajará con el time series conformada de ‘a_estacio’

## Series: a_estacio 
## ARIMA(2,0,1)(0,0,2)[12] with zero mean 
## 
## Coefficients:
##          ar1     ar2      ma1    sma1    sma2
##       0.4779  0.2408  -0.9531  0.2029  0.2801
## s.e.  0.0869  0.0855   0.0358  0.0825  0.1044
## 
## sigma^2 = 0.02311:  log likelihood = 79.27
## AIC=-146.54   AICc=-146.02   BIC=-127.8

En respuesta a la ARIMA(2,0,1)(0,0,2)[12]; se tiene que el modelo mide una parte no estacional con los componentes (2,0,1) y estacional (0,0,2)[12]; esto, permitirá capturar patrones tanto de corto plazo como de largo plazo.De esta manera, la primera parte contiene 2 componentes autoregresivos. 0 de estacionariedad y 1 con media móvil; por su parte, la segunda contiene una media móvil de 2 y el período de evaluación de 12 meses refiriendosea que los datos se presentan de manera anual.

Por su parte los valores bajos de AIC (Akaike Information Criterion) de -146.54 señala a complejidad del modelo, donde a menor valor mejor el modelo, AICc (Corrected AIC) de -146.02 y BIC (Bayesian Information Criterion) de -127.8 sugieren que el modelo se ajusta bien a los datos.

## numeric(0)

retorna numeric(0). Esto significa que no se detectaron puntos de cambio en la media de la serie temporal a_estacio utilizando el método ‘AMOC’. En otras palabras, la serie temporal no muestra evidencia de un cambio significativo en la media en ningún punto.

La salida numeric(0) indica que la media de la serie temporal a_estacio es constante a lo largo del tiempo, al menos según el método ‘AMOC’.

- Predicciones

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.85614, p-value = 1.487e-11

Se puede observar que que hay residuos extremos tanto en la parte inferior como superior de la linea de ajuste. Adicionalmente en el test de Shapiro el valor W: 0,856, establece que los datos se alejan de un comportamiento normal, en lina con este resultado el p-value: 1.487e-11 indica que existe suficiente evidencia estadistica para rechazar la hipotesis nula Ho:Los datos provienen de una distribucion normal.

## 
##  Box-Ljung test
## 
## data:  residuales
## X-squared = 0.001641, df = 1, p-value = 0.9677

Dado que la p-value es mayor a 0.05 se acepta la hipotesis nula, lo que significa que en primera instacia hay una buena especificación del modelo, las predicciones son confiables y los residuos son aleatorios, lo que permite al modelo capturar adecuadamente la estructura de dependencia temporal de los datos.

Modelo Prophet

Se realiza la aplicacion del modelo Prophet, debido a que el modelo solo reconoce las variables ds para unidad de tiempo y y para los valores de resultados se convierten las variables fechas y homicidios a este formato como se observa en la muestra a continuacion:

##           ds   y
## 1 2010-01-01 411
## 2 2010-02-01 392
## 3 2010-03-01 353
## 4 2010-04-01 389
## 5 2010-05-01 319
## 6 2010-06-01 397

Dado a que los datos se encuentran por mes, el ajuste del modelo se realiza activando la estacionalidad anual y desactivando la diaria y semanal.

El objetivo de prediccion para el modelo ajustado se establece para 11 periodos de la serie de tiempo, con una frecuencia mensual, lo que nos entregaria un resultado que abarca hasta diciembre 2024.

Predicciones

- Visualizacion del Modelo

Los puntos negros representan medidas reales de la serie de tiempo
La linea azul el pronóstico de Prophet el cual muestra un ajuste significativo a los datos reales
La banda azul representa el intervalo de incertidumbre, sin embargo este se muestra bastante conservador

En el grafico se puede observar el comportamiento durante el periodo de pandemia y la variación creciente del 2023 por fuera de las bandas predichas por el modelo, este comportamiento se muestra atipico en comparacion con el historico, lo que podria ser relevante analizar teniendo en cuenta la importancia de este indicador para la seguridad publica.

- Componentes del Modelo

A continuacion se presentan los componentes de modelo prophet en los que se observa la tendencia y los residuos, permitiendo identificar el comportamiento del modelo hasta la prediccion.

Para el caso de la tendencia, es claro que la prediccion presenta un comportamiento ascendente para el año 2024.

La grafica de residuales de la prediccion muestran, que a finales de septiembre los homicidios en el país tienden a la baja; adicionalmente, los meses de julio y enero son donde este indicador sufre un incremento significativo en comparacion con el resto de los meses.

Comparación entre modelos

Se crea una nueva serie de tiempo excluyendo un número de meses determinado con la finalidad de tomar ese tiempo como periodo de evaluación y comparación de los resultados de los modelos

- Creación del modelo Holt-Winter Train

Se realiza un nuevo modelo Holt-Winter para realizar el proceso de validacion de los modelos untilizando la serie de tiempo de entrenamiento con -6 meses.

## 
##  Augmented Dickey-Fuller Test
## 
## data:  nueva_ts
## Dickey-Fuller = 0.30191, Lag order = 5, p-value = 0.99
## alternative hypothesis: stationary

Dado que el p-value = 0.99 es mayor al nivel de significancia de 0.05 se rechaza la hipótesis alternativa de que la serie es estacionaria, luego de realizada la prueba Dicker - Fuller se identifica que se debe diferenciar la serie de tiempo 1 ves ya que no es estacionaria.

- Creación del modelo AutoArima Train

## Series: nueva_ts_estacionaria 
## ARIMA(0,0,1)(0,0,2)[12] with zero mean 
## 
## Coefficients:
##           ma1    sma1    sma2
##       -0.4869  0.1553  0.2768
## s.e.   0.0778  0.0803  0.1008
## 
## sigma^2 = 0.02123:  log likelihood = 82.49
## AIC=-156.97   AICc=-156.72   BIC=-144.62

En respuesta a la ARIMA(0,0,1)(0,0,2)[12]; se tiene que el modelo mide una parte no estacional con los componentes (0,0,1) y estacional (0,0,2)[12]; esto, permitirá capturar patrones tanto de corto plazo como de largo plazo. De esta manera, la primera parte contiene 0 componentes autoregresivos, 0 de estacionariedad y 1 con media móvil; por su parte, la segunda contiene una media móvil de 2 y el período de evaluación de 12 meses refiriendose a que los datos se presentan de manera anual.

- Creación del modelo Prophet Train

Las predicciones de los modelos ARIMA y Prophet se deben transformar para volver a la escala original, para esto se debe tener como referencia los datos sin transformar

- Valores Originales

Se imprime una muestra de los valores originales de la serie para los ultimos 6 meses, los cuales fueron utilizados para prueba del modelo

## actual_values
##  576  628  630 1177 1208 1264 
##    1    1    1    1    1    1

- Estimacion de Metricas de Modelo

## # A tibble: 3 × 9
##   Model           ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Holt-Winters  346.  482.  346.  56.2  56.2  2.09  1.56 0.478
## 2 ARIMA        -696.  760.  696. -73.1  73.1  4.20  2.45 0.444
## 3 Prophet       335.  500.  366.  56.0  58.6  2.21  1.61 0.484

## # A tibble: 3 × 9
##   Model           ME  RMSE   MAE   MPE  MAPE  MASE RMSSE  ACF1
##   <chr>        <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Holt-Winters  346.  482.  346. 56.2   56.2  2.09  1.56 0.478
## 2 ARIMA         252.  406.  318. 44.0   49.3  1.92  1.31 0.480
## 3 Prophet      -119.  349.  327. -1.27  35.4  1.97  1.13 0.463

Se realiza el ejercicio para diferentes periodos de evaluación: 12, 6 y 3 meses, es decir, se crea a partir de la serie de tiempo original una serie de entrenamiento y otra de prueba con el número de meses mencionado. En esta primera evaluación se usa la serie de datos trasnformados para que sea estacionaria antes de crear los modelos ARIMA y Prophet.

Para los 12 meses se obtuvieronlos siguientes resultados: